# install R packages

#install.packages("foreign")
#install.packages("randomizr")
#install.packages("readstata13")
#install.packages("AER")
#install.packages("nnet")
#install.packages("aod")
#install.packages("plyr")
#install.packages("descr")
#install.packages("xtable")
#install.packages("stargazer")
#install.packages("ri")
#install.packages("FSA")

setwd("") #Set working directory, in which Foos_John_PSRM_18_Jul_2016_load.dta was saved

library(foreign)  
library(randomizr)
library(readstata13)
library(AER) 
library(nnet) 
library(aod)     
library(plyr)
library(descr)
library(xtable)
library(stargazer)
library(ri)
library(FSA)
rm(list=ls())  

#Table 2

Raw.p<-c(.635, .008, .316) 

adjusted.p<-p.adjust(Raw.p, method="fdr")

data_upper<-as.data.frame(cbind(Raw.p, adjusted.p))

Raw.p<-c(.846, .582, .155, .130, .025, .018) 

adjusted.p<-p.adjust(Raw.p, method="fdr")

data_lower<-as.data.frame(cbind(Raw.p, adjusted.p))

Table_2<-rbind(data_upper, data_lower)

xtable(Table_2, digits=3)



#Table A1 

data<-read.dta13("Foos_John_PSRM_18_Jul_2016_load.dta")

sample<-subset(data, select = c(turnout2014, id, treat,  pid2,ward, pid, contact_hh, turnout2010, female, unknown, treat, gender))

sample$party_1[sample$pid2==1] <- 1
sample$party_1[sample$pid2==2] <- NA
sample$party_1[sample$pid2==3] <- NA
sample$party_1[sample$pid2==4] <- NA
sample$party_1[sample$pid2==5] <- NA
sample$party_1[sample$pid2==6] <- NA
sample$party_1[sample$pid2==7] <- NA

sample$party_2[sample$pid2==1] <- NA
sample$party_2[sample$pid2==2] <- 1
sample$party_2[sample$pid2==3] <- NA
sample$party_2[sample$pid2==4] <- NA
sample$party_2[sample$pid2==5] <- NA
sample$party_2[sample$pid2==6] <- NA
sample$party_2[sample$pid2==7] <- NA


sample$party_3[sample$pid2==1] <- NA
sample$party_3[sample$pid2==2] <- NA
sample$party_3[sample$pid2==3] <- 1
sample$party_3[sample$pid2==4] <- NA
sample$party_3[sample$pid2==5] <- NA
sample$party_3[sample$pid2==6] <- NA
sample$party_3[sample$pid2==7] <- NA

sample$party_4[sample$pid2==1] <- NA
sample$party_4[sample$pid2==2] <- NA
sample$party_4[sample$pid2==3] <- NA
sample$party_4[sample$pid2==4] <- 1
sample$party_4[sample$pid2==5] <- NA
sample$party_4[sample$pid2==6] <- NA
sample$party_4[sample$pid2==7] <- NA

sample$party_5[sample$pid2==1] <- NA
sample$party_5[sample$pid2==2] <- NA
sample$party_5[sample$pid2==3] <- NA
sample$party_5[sample$pid2==4] <- NA
sample$party_5[sample$pid2==5] <- 1
sample$party_5[sample$pid2==6] <- NA
sample$party_5[sample$pid2==7] <- NA

sample$party_6[sample$pid2==1] <- NA
sample$party_6[sample$pid2==2] <- NA
sample$party_6[sample$pid2==3] <- NA
sample$party_6[sample$pid2==4] <- NA
sample$party_6[sample$pid2==5] <- NA
sample$party_6[sample$pid2==6] <- 1
sample$party_6[sample$pid2==7] <- NA

sample$party_7[sample$pid2==1] <- NA
sample$party_7[sample$pid2==2] <- NA
sample$party_7[sample$pid2==3] <- NA
sample$party_7[sample$pid2==4] <- NA
sample$party_7[sample$pid2==5] <- NA
sample$party_7[sample$pid2==6] <- NA
sample$party_7[sample$pid2==7] <- 1


data$party_1[data$pid2==1] <- 1
data$party_1[data$pid2==2] <- 0
data$party_1[data$pid2==3] <- 0
data$party_1[data$pid2==4] <- 0
data$party_1[data$pid2==5] <- 0
data$party_1[data$pid2==6] <- 0
data$party_1[data$pid2==7] <- 0

data$party_2[data$pid2==1] <- 0
data$party_2[data$pid2==2] <- 1
data$party_2[data$pid2==3] <- 0
data$party_2[data$pid2==4] <- 0
data$party_2[data$pid2==5] <- 0
data$party_2[data$pid2==6] <- 0
data$party_2[data$pid2==7] <- 0

data$party_3[data$pid2==1] <- 0
data$party_3[data$pid2==2] <- 0
data$party_3[data$pid2==3] <- 1
data$party_3[data$pid2==4] <- 0
data$party_3[data$pid2==5] <- 0
data$party_3[data$pid2==6] <- 0
data$party_3[data$pid2==7] <- 0

data$party_4[data$pid2==1] <- 0
data$party_4[data$pid2==2] <- 0
data$party_4[data$pid2==3] <- 0
data$party_4[data$pid2==4] <- 1
data$party_4[data$pid2==5] <- 0
data$party_4[data$pid2==6] <- 0
data$party_4[data$pid2==7] <- 0

data$party_5[data$pid2==1] <- 0
data$party_5[data$pid2==2] <- 0
data$party_5[data$pid2==3] <- 0
data$party_5[data$pid2==4] <- 0
data$party_5[data$pid2==5] <- 1
data$party_5[data$pid2==6] <- 0
data$party_5[data$pid2==7] <- 0

data$party_6[data$pid2==1] <- 0
data$party_6[data$pid2==2] <- 0
data$party_6[data$pid2==3] <- 0
data$party_6[data$pid2==4] <- 0
data$party_6[data$pid2==5] <- 0
data$party_6[data$pid2==6] <- 1
data$party_6[data$pid2==7] <- 0

data$party_7[data$pid2==1] <- 0
data$party_7[data$pid2==2] <- 0
data$party_7[data$pid2==3] <- 0
data$party_7[data$pid2==4] <- 0
data$party_7[data$pid2==5] <- 0
data$party_7[data$pid2==6] <- 0
data$party_7[data$pid2==7] <- 1


sample$ward_1[sample$ward==0] <- 1
sample$ward_1[sample$ward==1] <- NA

sample$ward_2[sample$ward==0] <- NA
sample$ward_2[sample$ward==1] <- 1



library(descr)

sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_2))


#Balance Table 1 North

Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$contact_hh


N_control_North <-table(Z)[1]
N_canvass_North <-table(Z)[2]
N_leaflet_North <-table(Z)[3]

N_North=N_control_North+N_canvass_North+N_leaflet_North

X<-sample2$turnout2010

turnout_2010_North<-crosstab(Z, X, prop.r = TRUE)$prop.row

turnout_2010_North_control<-(round(turnout_2010_North[4], digits=3)*100)

turnout_2010_North_canvass<-(round(turnout_2010_North[5], digits=3)*100)

turnout_2010_North_leaflet<-(round(turnout_2010_North[6], digits=3)*100)


turnout_2010_North_overall<-cbind(turnout_2010_North_canvass, turnout_2010_North_leaflet, turnout_2010_North_control)

turnout_2010_North_overall


X<-sample2$female

female_North<-crosstab(Z, X, prop.r = TRUE)$prop.row

female_North_control<-(round(female_North[4], digits=3)*100)

female_North_canvass<-(round(female_North[5], digits=3)*100)

female_North_leaflet<-(round(female_North[6], digits=3)*100)

female_North_overall<-cbind(female_North_canvass, female_North_leaflet, female_North_control)

female_North_overall



X<-sample2$unknown

unknown_North<-crosstab(Z, X, prop.r = TRUE)$prop.row

unknown_North_control<-(round(unknown_North[4], digits=3)*100)

unknown_North_canvass<-(round(unknown_North[5], digits=3)*100)

unknown_North_leaflet<-(round(unknown_North[6], digits=3)*100)

unknown_North_overall<-cbind(unknown_North_canvass, unknown_North_leaflet, unknown_North_control)

unknown_North_overall


X<-sample2$pid2

pid_North<-crosstab(Z, X, prop.r = TRUE)$prop.row

Conservative_North_control<-(round(pid_North[1], digits=3)*100)

Conservative_North_canvass<-(round(pid_North[2], digits=3)*100)

Conservative_North_leaflet<-(round(pid_North[3], digits=3)*100)

AgainstCons_North_control<-(round(pid_North[4], digits=3)*100)

AgainstCons_North_canvass<-(round(pid_North[5], digits=3)*100)

AgainstCons_North_leaflet<-(round(pid_North[6], digits=3)*100)

Labour_North_control<-(round(pid_North[7], digits=3)*100)

Labour_North_canvass<-(round(pid_North[8], digits=3)*100)

Labour_North_leaflet<-(round(pid_North[9], digits=3)*100)

LibDem_North_control<-(round(pid_North[10], digits=3)*100)

LibDem_North_canvass<-(round(pid_North[11], digits=3)*100)

LibDem_North_leaflet<-(round(pid_North[12], digits=3)*100)

Undecided_North_control<-(round(pid_North[13], digits=3)*100)

Undecided_North_canvass<-(round(pid_North[14], digits=3)*100)

Undecided_North_leaflet<-(round(pid_North[15], digits=3)*100)

Nonvoter_North_control<-(round(pid_North[16], digits=3)*100)

Nonvoter_North_canvass<-(round(pid_North[17], digits=3)*100)

Nonvoter_North_leaflet<-(round(pid_North[18], digits=3)*100)

Missing_North_control<-(round(pid_North[19], digits=3)*100)

Missing_North_canvass<-(round(pid_North[20], digits=3)*100)

Missing_North_leaflet<-(round(pid_North[21], digits=3)*100)

Conservative_North<-cbind(Conservative_North_canvass, Conservative_North_leaflet, Conservative_North_control)

AgainstCons_North<-cbind(AgainstCons_North_canvass, AgainstCons_North_leaflet, AgainstCons_North_control)

Labour_North<-cbind(Labour_North_canvass, Labour_North_leaflet, Labour_North_control)

LibDem_North<-cbind(LibDem_North_canvass, LibDem_North_leaflet, LibDem_North_control)

Undecided_North<-cbind(Undecided_North_canvass, Undecided_North_leaflet, Undecided_North_control)

Nonvoter_North<-cbind(Nonvoter_North_canvass, Nonvoter_North_leaflet, Nonvoter_North_control)

Missing_North<-cbind(Missing_North_canvass, Missing_North_leaflet, Missing_North_control)


TableA1_North<-rbind(Conservative_North, AgainstCons_North, Labour_North, LibDem_North, Undecided_North, Nonvoter_North, Missing_North, female_North_overall, unknown_North_overall, turnout_2010_North_overall)



#Balance South

sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_1))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$compliance


N_control_South <-table(Z)[1]
N_canvass_South <-table(Z)[2]
N_leaflet_South <-table(Z)[3]

N_South=N_control_South+N_canvass_South+N_leaflet_South

N_South

X<-sample2$turnout2010

turnout_2010_South<-crosstab(Z, X, prop.r = TRUE)$prop.row

turnout_2010_South_control<-(round(turnout_2010_South[4], digits=3)*100)

turnout_2010_South_canvass<-(round(turnout_2010_South[5], digits=3)*100)

turnout_2010_South_leaflet<-(round(turnout_2010_South[6], digits=3)*100)


turnout_2010_South_overall<-cbind(turnout_2010_South_canvass, turnout_2010_South_leaflet, turnout_2010_South_control)

turnout_2010_South_overall


X<-sample2$female

female_South<-crosstab(Z, X, prop.r = TRUE)$prop.row

female_South_control<-(round(female_South[4], digits=3)*100)

female_South_canvass<-(round(female_South[5], digits=3)*100)

female_South_leaflet<-(round(female_South[6], digits=3)*100)

female_South_overall<-cbind(female_South_canvass, female_South_leaflet, female_South_control)

female_South_overall




X<-sample2$unknown

unknown_South<-crosstab(Z, X, prop.r = TRUE)$prop.row

unknown_South_control<-(round(unknown_South[4], digits=3)*100)

unknown_South_canvass<-(round(unknown_South[5], digits=3)*100)

unknown_South_leaflet<-(round(unknown_South[6], digits=3)*100)

unknown_South_overall<-cbind(unknown_South_canvass, unknown_South_leaflet, unknown_South_control)

unknown_South_overall


X<-sample2$pid2

pid_South<-crosstab(Z, X, prop.r = TRUE)$prop.row

Conservative_South_control<-(round(pid_South[1], digits=3)*100)

Conservative_South_canvass<-(round(pid_South[2], digits=3)*100)

Conservative_South_leaflet<-(round(pid_South[3], digits=3)*100)

AgainstCons_South_control<-(round(pid_South[4], digits=3)*100)

AgainstCons_South_canvass<-(round(pid_South[5], digits=3)*100)

AgainstCons_South_leaflet<-(round(pid_South[6], digits=3)*100)

Labour_South_control<-(round(pid_South[7], digits=3)*100)

Labour_South_canvass<-(round(pid_South[8], digits=3)*100)

Labour_South_leaflet<-(round(pid_South[9], digits=3)*100)

LibDem_South_control<-(round(pid_South[10], digits=3)*100)

LibDem_South_canvass<-(round(pid_South[11], digits=3)*100)

LibDem_South_leaflet<-(round(pid_South[12], digits=3)*100)

Undecided_South_control<-(round(pid_South[13], digits=3)*100)

Undecided_South_canvass<-(round(pid_South[14], digits=3)*100)

Undecided_South_leaflet<-(round(pid_South[15], digits=3)*100)

Nonvoter_South_control<-(round(pid_South[16], digits=3)*100)

Nonvoter_South_canvass<-(round(pid_South[17], digits=3)*100)

Nonvoter_South_leaflet<-(round(pid_South[18], digits=3)*100)

Missing_South_control<-(round(pid_South[19], digits=3)*100)

Missing_South_canvass<-(round(pid_South[20], digits=3)*100)

Missing_South_leaflet<-(round(pid_South[21], digits=3)*100)

Conservative_South<-cbind(Conservative_South_canvass, Conservative_South_leaflet, Conservative_South_control)

AgainstCons_South<-cbind(AgainstCons_South_canvass, AgainstCons_South_leaflet, AgainstCons_South_control)

Labour_South<-cbind(Labour_South_canvass, Labour_South_leaflet, Labour_South_control)

LibDem_South<-cbind(LibDem_South_canvass, LibDem_South_leaflet, LibDem_South_control)

Undecided_South<-cbind(Undecided_South_canvass, Undecided_South_leaflet, Undecided_South_control)

Nonvoter_South<-cbind(Nonvoter_South_canvass, Nonvoter_South_leaflet, Nonvoter_South_control)

Missing_South<-cbind(Missing_South_canvass, Missing_South_leaflet, Missing_South_control)




TableA1_South<-rbind(Conservative_South, AgainstCons_South, Labour_South, LibDem_South, Undecided_South, Nonvoter_South, Missing_South, female_South_overall, unknown_South_overall, turnout_2010_South_overall)

TableA1_South


TableA1<-rbind(TableA1_North, N_North, TableA1_South, N_South)

TableA1

colnames(TableA1)[1] <- "Canvass"
colnames(TableA1)[2] <- "Leaflet"
colnames(TableA1)[3] <- "Control"
rownames(TableA1)[1] <- "Conservative - North"
rownames(TableA1)[2] <- "Against Conservatives- North"
rownames(TableA1)[3] <- "Labour - North"
rownames(TableA1)[4] <- "LibDem - North"
rownames(TableA1)[5] <- "Undecided - North"
rownames(TableA1)[6] <- "Non-voter - North"
rownames(TableA1)[7] <- "Missing - North"
rownames(TableA1)[8] <- "Gender Female - North"
rownames(TableA1)[9] <- "Gender Unknown - North"
rownames(TableA1)[10] <- "Turnout 2010 - North"
rownames(TableA1)[11] <- "N North"
rownames(TableA1)[12] <- "Conservative - South"
rownames(TableA1)[13] <- "Against Conservatives- South"
rownames(TableA1)[14] <- "Labour - South"
rownames(TableA1)[15] <- "LibDem - South"
rownames(TableA1)[16] <- "Undecided - South"
rownames(TableA1)[17] <- "Non-voter - South"
rownames(TableA1)[18] <- "Missing - South"
rownames(TableA1)[19] <- "Gender Female - South"
rownames(TableA1)[20] <- "Gender Unknown - South"
rownames(TableA1)[21] <- "Turnout 2010 - South"
rownames(TableA1)[22] <- "N South"


digits<-matrix(1,nrow=22,ncol=4)

digits[c(11,22),]<-0

xtable(TableA1, digits = digits)



#Table A2 

#Turnout North

sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_2))

Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$contact_hh

turnout_2014_North<-crosstab(Z, Y, prop.r = TRUE)$prop.row

turnout_2014_North_control<-(round(turnout_2014_North[4], digits=3)*100)

turnout_2014_North_canvass<-(round(turnout_2014_North[5], digits=3)*100)

turnout_2014_North_leaflet<-(round(turnout_2014_North[6], digits=3)*100)

turnout_2014_North_combined<-(round(prop.table(table(Y))[2], digits=3)*100)

turnout_2014_North_overall<-cbind(turnout_2014_North_canvass, turnout_2014_North_leaflet, turnout_2014_North_control, turnout_2014_North_combined)

turnout_2014_North_overall


#N North combined


N_control_North<-table(Z)[1]

N_canvass_North<-table(Z)[2]

N_leaflet_North<-table(Z)[3]

N_combined_North<-(margin.table(table(Z)))

N_combined_North_combined<-cbind(N_canvass_North, N_leaflet_North, N_control_North, N_combined_North)


sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_2))

sample2<-subset(sample2,!is.na(party_1))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$contact_hh


turnout_2014_Cons_North<-crosstab(Z, Y, prop.r = TRUE)$prop.row

turnout_2014_Cons_North_control<-(round(turnout_2014_Cons_North[4], digits=3)*100)

turnout_2014_Cons_North_canvass<-(round(turnout_2014_Cons_North[5], digits=3)*100)

turnout_2014_Cons_North_leaflet<-(round(turnout_2014_Cons_North[6], digits=3)*100)

turnout_2014_Cons_North_combined<-(round(prop.table(table(Y))[2], digits=3)*100)


turnout_2014_Cons_North_overall<-cbind(turnout_2014_Cons_North_canvass, turnout_2014_Cons_North_leaflet, turnout_2014_Cons_North_control, turnout_2014_Cons_North_combined)

turnout_2014_Cons_North_overall



#N Cons North

N_Cons_control_North<-table(Z)[1]

N_Cons_canvass_North<-table(Z)[2]

N_Cons_leaflet_North<-table(Z)[3]

N_Cons_combined_North<-(margin.table(table(Z)))

N_Cons_North<-cbind(N_Cons_canvass_North, N_Cons_leaflet_North, N_Cons_control_North, N_Cons_combined_North)




sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_2))

sample2<-subset(sample2,!is.na(party_2))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$contact_hh

turnout_2014_AgainstCons_North<-crosstab(Z, Y, prop.r = TRUE)$prop.row


turnout_2014_AgainstCons_North_control<-(round(turnout_2014_AgainstCons_North[4], digits=3)*100)

turnout_2014_AgainstCons_North_canvass<-(round(turnout_2014_AgainstCons_North[5], digits=3)*100)

turnout_2014_AgainstCons_North_leaflet<-(round(turnout_2014_AgainstCons_North[6], digits=3)*100)

turnout_2014_AgainstCons_North_combined<-(round(prop.table(table(Y))[2], digits=3)*100)


turnout_2014_AgainstCons_North_overall<-cbind(turnout_2014_AgainstCons_North_canvass, turnout_2014_AgainstCons_North_leaflet, turnout_2014_AgainstCons_North_control, turnout_2014_AgainstCons_North_combined)

turnout_2014_AgainstCons_North_overall


#N Against Cons North

N_AgainstCons_control_North<-table(Z)[1]

N_AgainstCons_canvass_North<-table(Z)[2]

N_AgainstCons_leaflet_North<-table(Z)[3]

N_AgainstCons_combined_North<-(margin.table(table(Z)))

N_AgainstCons_North<-cbind(N_AgainstCons_canvass_North, N_AgainstCons_leaflet_North, N_AgainstCons_control_North, N_AgainstCons_combined_North)




sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_2))

sample2<-subset(sample2,!is.na(party_3))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$contact_hh

turnout_2014_Labour_North<-crosstab(Z, Y, prop.r = TRUE)$prop.row


turnout_2014_Labour_North_control<-(round(turnout_2014_Labour_North[4], digits=3)*100)

turnout_2014_Labour_North_canvass<-(round(turnout_2014_Labour_North[5], digits=3)*100)

turnout_2014_Labour_North_leaflet<-(round(turnout_2014_Labour_North[6], digits=3)*100)

turnout_2014_Labour_North_combined<-(round(prop.table(table(Y))[2], digits=3)*100)



turnout_2014_Labour_North_overall<-cbind(turnout_2014_Labour_North_canvass, turnout_2014_Labour_North_leaflet, turnout_2014_Labour_North_control, turnout_2014_Labour_North_combined)

turnout_2014_Labour_North_overall


#N Labour North

N_Labour_control_North<-table(Z)[1]

N_Labour_canvass_North<-table(Z)[2]

N_Labour_leaflet_North<-table(Z)[3]

N_Labour_combined_North<-(margin.table(table(Z)))

N_Labour_North<-cbind(N_Labour_canvass_North, N_Labour_leaflet_North, N_Labour_control_North, N_Labour_combined_North)





sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_2))

sample2<-subset(sample2,!is.na(party_4))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$contact_hh


turnout_2014_LibDem_North<-crosstab(Z, Y, prop.r = TRUE)$prop.row


turnout_2014_LibDem_North_control<-(round(turnout_2014_LibDem_North[4], digits=3)*100)

turnout_2014_LibDem_North_canvass<-(round(turnout_2014_LibDem_North[5], digits=3)*100)

turnout_2014_LibDem_North_leaflet<-(round(turnout_2014_LibDem_North[6], digits=3)*100)

turnout_2014_LibDem_North_combined<-(round(prop.table(table(Y))[2], digits=3)*100)


turnout_2014_LibDem_North_overall<-cbind(turnout_2014_LibDem_North_canvass, turnout_2014_LibDem_North_leaflet, turnout_2014_LibDem_North_control, turnout_2014_LibDem_North_combined)

turnout_2014_LibDem_North_overall


#N LibDem North

N_LibDem_control_North<-table(Z)[1]

N_LibDem_canvass_North<-table(Z)[2]

N_LibDem_leaflet_North<-table(Z)[3]

N_LibDem_combined_North<-(margin.table(table(Z)))

N_LibDem_North<-cbind(N_LibDem_canvass_North, N_LibDem_leaflet_North, N_LibDem_control_North, N_LibDem_combined_North)




sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_2))

sample2<-subset(sample2,!is.na(party_5))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$contact_hh


turnout_2014_Undecided_North<-crosstab(Z, Y, prop.r = TRUE)$prop.row

turnout_2014_Undecided_North_control<-(round(turnout_2014_Undecided_North[4], digits=3)*100)

turnout_2014_Undecided_North_canvass<-(round(turnout_2014_Undecided_North[5], digits=3)*100)

turnout_2014_Undecided_North_leaflet<-(round(turnout_2014_Undecided_North[6], digits=3)*100)

turnout_2014_Undecided_North_combined<-(round(prop.table(table(Y))[2], digits=3)*100)


turnout_2014_Undecided_North_overall<-cbind(turnout_2014_Undecided_North_canvass, turnout_2014_Undecided_North_leaflet, turnout_2014_Undecided_North_control, turnout_2014_Undecided_North_combined)

turnout_2014_Undecided_North_overall


#N Undecided North

N_Undecided_control_North<-table(Z)[1]

N_Undecided_canvass_North<-table(Z)[2]

N_Undecided_leaflet_North<-table(Z)[3]

N_Undecided_combined_North<-(margin.table(table(Z)))

N_Undecided_North<-cbind(N_Undecided_canvass_North, N_Undecided_leaflet_North, N_Undecided_control_North, N_Undecided_combined_North)



sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_2))

sample2<-subset(sample2,!is.na(party_6))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$contact_hh


turnout_2014_Nonvoter_North<-crosstab(Z, Y, prop.r = TRUE)$prop.row

turnout_2014_Nonvoter_North_control<-(round(turnout_2014_Nonvoter_North[4], digits=3)*100)

turnout_2014_Nonvoter_North_canvass<-(round(turnout_2014_Nonvoter_North[5], digits=3)*100)

turnout_2014_Nonvoter_North_leaflet<-(round(turnout_2014_Nonvoter_North[6], digits=3)*100)

turnout_2014_Nonvoter_North_combined<-(round(prop.table(table(Y))[2], digits=3)*100)


turnout_2014_Nonvoter_North_overall<-cbind(turnout_2014_Nonvoter_North_canvass, turnout_2014_Nonvoter_North_leaflet, turnout_2014_Nonvoter_North_control, turnout_2014_Nonvoter_North_combined)

turnout_2014_Nonvoter_North_overall



#N Non-voter North

N_Nonvoter_control_North<-table(Z)[1]

N_Nonvoter_canvass_North<-table(Z)[2]

N_Nonvoter_leaflet_North<-table(Z)[3]

N_Nonvoter_combined_North<-(margin.table(table(Z)))

N_Nonvoter_North<-cbind(N_Nonvoter_canvass_North, N_Nonvoter_leaflet_North, N_Nonvoter_control_North, N_Nonvoter_combined_North)




sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_2))

sample2<-subset(sample2,!is.na(party_7))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$contact_hh


turnout_2014_Missing_North<-crosstab(Z, Y, prop.r = TRUE)$prop.row

turnout_2014_Missing_North_control<-(round(turnout_2014_Missing_North[4], digits=3)*100)

turnout_2014_Missing_North_canvass<-(round(turnout_2014_Missing_North[5], digits=3)*100)

turnout_2014_Missing_North_leaflet<-(round(turnout_2014_Missing_North[6], digits=3)*100)

turnout_2014_Missing_North_combined<-(round(prop.table(table(Y))[2], digits=3)*100)


turnout_2014_Missing_North_overall<-cbind(turnout_2014_Missing_North_canvass, turnout_2014_Missing_North_leaflet, turnout_2014_Missing_North_control, turnout_2014_Missing_North_combined)

turnout_2014_Missing_North_overall



#N Missing North

N_Missing_control_North<-table(Z)[1]

N_Missing_canvass_North<-table(Z)[2]

N_Missing_leaflet_North<-table(Z)[3]

N_Missing_combined_North<-(margin.table(table(Z)))

N_Missing_North<-cbind(N_Missing_canvass_North, N_Missing_leaflet_North, N_Missing_control_North, N_Missing_combined_North)




#Contact

sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_2))

Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$contact_hh


contact_North<-crosstab(Z, D, prop.r = TRUE)$prop.row


contact_North_control<-(round(contact_North[4], digits=3)*100)

contact_North_canvass<-(round(contact_North[5], digits=3)*100)

contact_North_leaflet<-(round(contact_North[6], digits=3)*100)

turnout_contact_North_combined<-(round(prop.table(table(D))[2], digits=3)*100)



contact_North_overall<-cbind(contact_North_canvass, contact_North_leaflet, contact_North_control, turnout_contact_North_combined)

contact_North_overall



TableA2_North<-rbind(turnout_2014_North_overall, N_combined_North_combined, turnout_2014_Cons_North_overall, N_Cons_North, turnout_2014_AgainstCons_North_overall, N_AgainstCons_North, turnout_2014_Labour_North_overall, N_Labour_North, turnout_2014_LibDem_North_overall, N_LibDem_North, turnout_2014_Undecided_North_overall, N_Undecided_North, turnout_2014_Nonvoter_North_overall, N_Nonvoter_North, turnout_2014_Missing_North_overall, N_Missing_North, contact_North_overall)

TableA2_North



#Turnout South


#####

sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_1))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$compliance


turnout_2014_South<-crosstab(Z, Y, prop.r = TRUE)$prop.row

turnout_2014_South_control<-(round(turnout_2014_South[4], digits=3)*100)

turnout_2014_South_canvass<-(round(turnout_2014_South[5], digits=3)*100)

turnout_2014_South_leaflet<-(round(turnout_2014_South[6], digits=3)*100)

turnout_2014_South_combined<-(round(prop.table(table(Y))[2], digits=3)*100)


turnout_2014_South_overall<-cbind(turnout_2014_South_canvass, turnout_2014_South_leaflet, turnout_2014_South_control, turnout_2014_South_combined)

turnout_2014_South_overall


#N South Combined

N_combined_control_South<-table(Z)[1]

N_combined_canvass_South<-table(Z)[2]

N_combined_leaflet_South<-table(Z)[3]

N_combined_South<-(margin.table(table(Z)))

N_combined_South_combined<-cbind(N_combined_canvass_South, N_combined_leaflet_South, N_combined_control_South, N_combined_South)





sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_1))

sample2<-subset(sample2,!is.na(party_1))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$compliance


turnout_2014_Cons_South<-crosstab(Z, Y, prop.r = TRUE)$prop.row


turnout_2014_Cons_South_control<-(round(turnout_2014_Cons_South[4], digits=3)*100)

turnout_2014_Cons_South_canvass<-(round(turnout_2014_Cons_South[5], digits=3)*100)

turnout_2014_Cons_South_leaflet<-(round(turnout_2014_Cons_South[6], digits=3)*100)

turnout_2014_Cons_South_combined<-(round(prop.table(table(Y))[2], digits=3)*100)



turnout_2014_Cons_South_overall<-cbind(turnout_2014_Cons_South_canvass, turnout_2014_Cons_South_leaflet, turnout_2014_Cons_South_control, turnout_2014_Cons_South_combined)

turnout_2014_Cons_South_overall


#N South Cons

N_Cons_control_South<-table(Z)[1]

N_Cons_canvass_South<-table(Z)[2]

N_Cons_leaflet_South<-table(Z)[3]

N_Cons_combined_South<-(margin.table(table(Z)))

N_Cons_South<-cbind(N_Cons_canvass_South, N_Cons_leaflet_South, N_Cons_control_South, N_Cons_combined_South)



sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_1))

sample2<-subset(sample2,!is.na(party_2))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$compliance


turnout_2014_AgainstCons_South<-crosstab(Z, Y, prop.r = TRUE)$prop.row


turnout_2014_AgainstCons_South_control<-(round(turnout_2014_AgainstCons_South[4], digits=3)*100)

turnout_2014_AgainstCons_South_canvass<-(round(turnout_2014_AgainstCons_South[5], digits=3)*100)

turnout_2014_AgainstCons_South_leaflet<-(round(turnout_2014_AgainstCons_South[6], digits=3)*100)

turnout_2014_AgainstCons_South_combined<-(round(prop.table(table(Y))[2], digits=3)*100)


turnout_2014_AgainstCons_South_overall<-cbind(turnout_2014_AgainstCons_South_canvass, turnout_2014_AgainstCons_South_leaflet, turnout_2014_AgainstCons_South_control, turnout_2014_AgainstCons_South_combined)

turnout_2014_AgainstCons_South_overall




#N South Against Cons

N_AgainstCons_control_South<-table(Z)[1]

N_AgainstCons_canvass_South<-table(Z)[2]

N_AgainstCons_leaflet_South<-table(Z)[3]

N_AgainstCons_combined_South<-(margin.table(table(Z)))

N_AgainstCons_South<-cbind(N_AgainstCons_canvass_South, N_AgainstCons_leaflet_South, N_AgainstCons_control_South, N_AgainstCons_combined_South)




sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_1))

sample2<-subset(sample2,!is.na(party_3))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$compliance


turnout_2014_Labour_South<-crosstab(Z, Y, prop.r = TRUE)$prop.row


turnout_2014_Labour_South_control<-(round(turnout_2014_Labour_South[4], digits=3)*100)

turnout_2014_Labour_South_canvass<-(round(turnout_2014_Labour_South[5], digits=3)*100)

turnout_2014_Labour_South_leaflet<-(round(turnout_2014_Labour_South[6], digits=3)*100)

turnout_2014_Labour_South_combined<-(round(prop.table(table(Y))[2], digits=3)*100)



turnout_2014_Labour_South_overall<-cbind(turnout_2014_Labour_South_canvass, turnout_2014_Labour_South_leaflet, turnout_2014_Labour_South_control, turnout_2014_Labour_South_combined)

turnout_2014_Labour_South_overall




#N South Labour

N_Labour_control_South<-table(Z)[1]

N_Labour_canvass_South<-table(Z)[2]

N_Labour_leaflet_South<-table(Z)[3]

N_Labour_combined_South<-(margin.table(table(Z)))

N_Labour_South<-cbind(N_Labour_canvass_South, N_Labour_leaflet_South, N_Labour_control_South, N_Labour_combined_South)





sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_1))

sample2<-subset(sample2,!is.na(party_4))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$compliance


turnout_2014_LibDem_South<-crosstab(Z, Y, prop.r = TRUE)$prop.row


turnout_2014_LibDem_South_control<-(round(turnout_2014_LibDem_South[4], digits=3)*100)

turnout_2014_LibDem_South_canvass<-(round(turnout_2014_LibDem_South[5], digits=3)*100)

turnout_2014_LibDem_South_leaflet<-(round(turnout_2014_LibDem_South[6], digits=3)*100)

turnout_2014_LibDem_South_combined<-(round(prop.table(table(Y))[2], digits=3)*100)


turnout_2014_LibDem_South_overall<-cbind(turnout_2014_LibDem_South_canvass, turnout_2014_LibDem_South_leaflet, turnout_2014_LibDem_South_control, turnout_2014_LibDem_South_combined)

turnout_2014_LibDem_South_overall




#N South LibDem

N_LibDem_control_South<-table(Z)[1]

N_LibDem_canvass_South<-table(Z)[2]

N_LibDem_leaflet_South<-table(Z)[3]

N_LibDem_combined_South<-(margin.table(table(Z)))

N_LibDem_South<-cbind(N_LibDem_canvass_South, N_LibDem_leaflet_South, N_LibDem_control_South, N_LibDem_combined_South)




sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_1))

sample2<-subset(sample2,!is.na(party_5))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$compliance



turnout_2014_Undecided_South<-crosstab(Z, Y, prop.r = TRUE)$prop.row


turnout_2014_Undecided_South_control<-(round(turnout_2014_Undecided_South[4], digits=3)*100)

turnout_2014_Undecided_South_canvass<-(round(turnout_2014_Undecided_South[5], digits=3)*100)

turnout_2014_Undecided_South_leaflet<-(round(turnout_2014_Undecided_South[6], digits=3)*100)

turnout_2014_Undecided_South_combined<-(round(prop.table(table(Y))[2], digits=3)*100)



turnout_2014_Undecided_South_overall<-cbind(turnout_2014_Undecided_South_canvass, turnout_2014_Undecided_South_leaflet, turnout_2014_Undecided_South_control, turnout_2014_Undecided_South_combined)

turnout_2014_Undecided_South_overall



#N South Undecided

N_Undecided_control_South<-table(Z)[1]

N_Undecided_canvass_South<-table(Z)[2]

N_Undecided_leaflet_South<-table(Z)[3]

N_Undecided_combined_South<-(margin.table(table(Z)))

N_Undecided_South<-cbind(N_Undecided_canvass_South, N_Undecided_leaflet_South, N_Undecided_control_South, N_Undecided_combined_South)




sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_1))

sample2<-subset(sample2,!is.na(party_6))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$compliance




turnout_2014_Nonvoter_South<-crosstab(Z, Y, prop.r = TRUE)$prop.row


turnout_2014_Nonvoter_South_control<-(round(turnout_2014_Nonvoter_South[4], digits=3)*100)

turnout_2014_Nonvoter_South_canvass<-(round(turnout_2014_Nonvoter_South[5], digits=3)*100)

turnout_2014_Nonvoter_South_leaflet<-(round(turnout_2014_Nonvoter_South[6], digits=3)*100)

turnout_2014_Nonvoter_South_combined<-(round(prop.table(table(Y))[2], digits=3)*100)


turnout_2014_Nonvoter_South_overall<-cbind(turnout_2014_Nonvoter_South_canvass, turnout_2014_Nonvoter_South_leaflet, turnout_2014_Nonvoter_South_control, turnout_2014_Nonvoter_South_combined)

turnout_2014_Nonvoter_South_overall


#N South Non-voter

N_Nonvoter_control_South<-table(Z)[1]

N_Nonvoter_canvass_South<-table(Z)[2]

N_Nonvoter_leaflet_South<-table(Z)[3]

N_Nonvoter_combined_South<-(margin.table(table(Z)))

N_Nonvoter_South<-cbind(N_Nonvoter_canvass_South, N_Nonvoter_leaflet_South, N_Nonvoter_control_South, N_Nonvoter_combined_South)




sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_1))

sample2<-subset(sample2,!is.na(party_7))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$compliance




turnout_2014_Missing_South<-crosstab(Z, Y, prop.r = TRUE)$prop.row


turnout_2014_Missing_South_control<-(round(turnout_2014_Missing_South[4], digits=3)*100)

turnout_2014_Missing_South_canvass<-(round(turnout_2014_Missing_South[5], digits=3)*100)

turnout_2014_Missing_South_leaflet<-(round(turnout_2014_Missing_South[6], digits=3)*100)

turnout_2014_Missing_South_combined<-(round(prop.table(table(Y))[2], digits=3)*100)


turnout_2014_Missing_South_overall<-cbind(turnout_2014_Missing_South_canvass, turnout_2014_Missing_South_leaflet, turnout_2014_Missing_South_control, turnout_2014_Missing_South_combined)

turnout_2014_Missing_South_overall


#N South Missing

N_Missing_control_South<-table(Z)[1]

N_Missing_canvass_South<-table(Z)[2]

N_Missing_leaflet_South<-table(Z)[3]

N_Missing_combined_South<-(margin.table(table(Z)))

N_Missing_South<-cbind(N_Missing_canvass_South, N_Missing_leaflet_South, N_Missing_control_South, N_Missing_combined_South)




sample<-subset(sample,!is.na(id))

sample1<-subset(sample,!is.na(treat))

sample2<-subset(sample1,!is.na(ward_1))


Y<-sample2$turnout2014
Z<-sample2$treat
D<-sample2$contact_hh

contact_South<-crosstab(Z, D, prop.r = TRUE)$prop.row


contact_South_control<-(round(contact_South[4], digits=3)*100)

contact_South_canvass<-(round(contact_South[5], digits=3)*100)

contact_South_leaflet<-(round(contact_South[6], digits=3)*100)

contact_South_combined<-(round(prop.table(table(D))[2], digits=3)*100)


contact_South_overall<-cbind(contact_South_canvass, contact_South_leaflet, contact_South_control, contact_South_combined)

contact_South_overall



TableA2_South<-rbind(turnout_2014_South_overall,  N_combined_South, turnout_2014_Cons_South_overall, N_Cons_South, turnout_2014_AgainstCons_South_overall, N_AgainstCons_South, turnout_2014_Labour_South_overall, N_Labour_South, turnout_2014_LibDem_South_overall, N_LibDem_South, turnout_2014_Undecided_South_overall, N_Undecided_South, turnout_2014_Nonvoter_South_overall, N_Nonvoter_South, turnout_2014_Missing_South_overall, N_Missing_South, contact_South_overall)

TableA2_South


TableA2<-rbind(TableA2_North, TableA2_South)

TableA2

colnames(TableA2)[1] <- "Canvass"
colnames(TableA2)[2] <- "Leaflet"
colnames(TableA2)[3] <- "Control"
colnames(TableA2)[4] <- "Combined"
rownames(TableA2)[1] <- "Combined - North"
rownames(TableA2)[2] <- "N Combined - North"
rownames(TableA2)[3] <- "Conservative - North"
rownames(TableA2)[4] <- "N Conservative North"
rownames(TableA2)[5] <- "Against Conservatives- North"
rownames(TableA2)[6] <- "N Against Conservatives- North"
rownames(TableA2)[7] <- "Labour - North"
rownames(TableA2)[8] <- "N Labour - North"
rownames(TableA2)[9] <- "LibDem - North"
rownames(TableA2)[10] <- "N LibDem - North"
rownames(TableA2)[11] <- "Undecided - North"
rownames(TableA2)[12] <- "N Undecided - North"
rownames(TableA2)[13] <- "Non-voter - North"
rownames(TableA2)[14] <- "N Non-voter - North"
rownames(TableA2)[15] <- "Missing - North"
rownames(TableA2)[16] <- "N Missing - North"
rownames(TableA2)[17] <- "Contact - North"
rownames(TableA2)[18] <- "Combined - South"
rownames(TableA2)[19] <- "N Combined - South"
rownames(TableA2)[20] <- "Conservative - South"
rownames(TableA2)[21] <- "N Conservative South"
rownames(TableA2)[22] <- "Against Conservative - South"
rownames(TableA2)[23] <- "N Against Conservative South"
rownames(TableA2)[24] <- "Labour - South"
rownames(TableA2)[25] <- "N Labour - South"
rownames(TableA2)[26] <- "LibDem - South"
rownames(TableA2)[27] <- "N LibDem - South"
rownames(TableA2)[28] <- "Undecided - South"
rownames(TableA2)[29] <- "N Undecided - South"
rownames(TableA2)[30] <- "Non-voter - South"
rownames(TableA2)[31] <- "N Non-voter - South"
rownames(TableA2)[32] <- "Missing - South"
rownames(TableA2)[33] <- "N Missing - South"
rownames(TableA2)[34] <- "Contact - South"


digits<-matrix(1,nrow=34,ncol=5)

digits[c(2,4,6,8,10,12,14,16,19,21,23,25,27,29,31,33),]<-0

xtable(TableA2, digits = digits)



#Figure A1

data<-read.dta13("Foos_John_PSRM_18_Jul_2016_load.dta")

data$party_1[data$pid2==1] <- 1
data$party_1[data$pid2==2] <- 0
data$party_1[data$pid2==3] <- 0
data$party_1[data$pid2==4] <- 0
data$party_1[data$pid2==5] <- 0
data$party_1[data$pid2==6] <- 0
data$party_1[data$pid2==7] <- 0


data$party_2[data$pid2==1] <- 0
data$party_2[data$pid2==2] <- 1
data$party_2[data$pid2==3] <- 0
data$party_2[data$pid2==4] <- 0
data$party_2[data$pid2==5] <- 0
data$party_2[data$pid2==6] <- 0
data$party_2[data$pid2==7] <- 0

data$party_3[data$pid2==1] <- 0
data$party_3[data$pid2==2] <- 0
data$party_3[data$pid2==3] <- 1
data$party_3[data$pid2==4] <- 0
data$party_3[data$pid2==5] <- 0
data$party_3[data$pid2==6] <- 0
data$party_3[data$pid2==7] <- 0

data$party_4[data$pid2==1] <- 0
data$party_4[data$pid2==2] <- 0
data$party_4[data$pid2==3] <- 0
data$party_4[data$pid2==4] <- 1
data$party_4[data$pid2==5] <- 0
data$party_4[data$pid2==6] <- 0
data$party_4[data$pid2==7] <- 0

data$party_5[data$pid2==1] <- 0
data$party_5[data$pid2==2] <- 0
data$party_5[data$pid2==3] <- 0
data$party_5[data$pid2==4] <- 0
data$party_5[data$pid2==5] <- 1
data$party_5[data$pid2==6] <- 0
data$party_5[data$pid2==7] <- 0

data$party_6[data$pid2==1] <- 0
data$party_6[data$pid2==2] <- 0
data$party_6[data$pid2==3] <- 0
data$party_6[data$pid2==4] <- 0
data$party_6[data$pid2==5] <- 0
data$party_6[data$pid2==6] <- 1
data$party_6[data$pid2==7] <- 0

data$party_7[data$pid2==1] <- 0
data$party_7[data$pid2==2] <- 0
data$party_7[data$pid2==3] <- 0
data$party_7[data$pid2==4] <- 0
data$party_7[data$pid2==5] <- 0
data$party_7[data$pid2==6] <- 0
data$party_7[data$pid2==7] <- 1

data$treat[data$group==0] <- 1
data$treat[data$group==1] <- 2
data$treat[data$group==2] <- 3


#Get household compliance

data1<-subset(data, select = c(id, treat, ward, turnout2010, female, unknown, party_2, party_3, party_4, party_5, party_6, party_7, party_code, contact_hh))

individual_level <- data1
household_level <- ddply(individual_level, c("id", "ward"), summarize, treatment=mean(contact_hh))   
with(household_level, table(ward, treatment))  ## This is how many households are assigned to treatment.

#Balance

#Execute code and save p-value and figure

data1<-subset(data, select = c(id, treat, ward, turnout2010, female, unknown, party_2, party_3, party_4, party_5, party_6, party_7))


individual_level <- data1
household_level <- ddply(individual_level, c("id", "ward"), summarize, treatment=mean(treat))   
with(household_level, table(ward, treatment))  ## This is how many treatments to allocate.

Somerset_ra <- function(){
  household_level <- within(household_level,{
    Z_sim <- rep(NA, nrow(household_level))
    Z_sim[ward==0] <- complete_ra(N = sum(ward==0), m_each=c(559, 500, 500))
    Z_sim[ward==1] <- complete_ra(N = sum(ward==1), m_each=c(660, 500, 500))
  })
  individual_level <- join(individual_level, household_level, by = c("ward", "id"))
  return(individual_level$Z_sim)
}

table(Somerset_ra())


set.seed(1234567) 
nsims <- 5000

Z<-individual_level$treat #Three-armed experiment (1 = control, 2 = treatment group 1, 3= treatment group 2)
blockvar<-individual_level$ward # blocks = 2 experimental sides, 1 = "South", 0 = "North"
clustvar<-individual_level$id # cluster random assignment at the household level, id=household identifier
X <- as.matrix(individual_level[,4:12])

loglikstore <- rep(NA,nsims)

for (i in 1:nsims) {
  loglikstore[i] <- summary(multinom(Somerset_ra()~X))$deviance/2 
}

loglik <- summary(multinom(Z~X))$deviance/2

mean(loglik <= loglikstore) # 0.7924

hist(loglikstore)
abline(v=loglik, col="red", lwd=2)

pdf(paste("balance.pdf"),w=6,h=5)
hist(loglikstore, breaks = 100, main = paste("Sampling distribution of simulated log likelihoods"), xlab=("Log Likelihood"))
abline(v=loglik, col="red", lwd=3)
dev.off() 


#########################
#Interaction balance

######################

data1<-subset(data, select = c(id, treat, ward, turnout2010, female, unknown, pid))

individual_level <- data1
## better not to have r objects that are named the same as functions.
household_level <- ddply(individual_level, c("id", "ward"), summarize, treatment=mean(treat))   
with(household_level, table(ward, treatment))  ## This is how many treatments to allocate.

Somerset_ra <- function(){
  household_level <- within(household_level,{
    Z_sim <- rep(NA, nrow(household_level))
    Z_sim[ward==0] <- complete_ra(N = sum(ward==0), m_each=c(559, 500, 500))
    Z_sim[ward==1] <- complete_ra(N = sum(ward==1), m_each=c(660, 500, 500))
  })
  individual_level <- join(individual_level, household_level, by = c("ward", "id"))
  return(individual_level$Z_sim)
}

table(Somerset_ra())

set.seed(1234567) 
nsims <- 5000

Z<-individual_level$treat #Three-armed experiment (1 = control, 2 = treatment group 1, 3= treatment group 2)
blockvar<-individual_level$block # blocks = 2 experimental sides, 1 = "South", 0 = "North"
clustvar<-individual_level$id # cluster random assignment at the household level, id=household identifier
X <- as.matrix(individual_level[,4:6])
pid<-individual_level$pid

loglikstore <- rep(NA,nsims)

for (i in 1:nsims) {
  loglikstore[i] <- summary(multinom(Somerset_ra()~factor(pid) + X + factor(pid)*X))$deviance/2 
}

loglik <- summary(multinom(Z~factor(pid) + X + factor(pid)*X))$deviance/2

mean(loglik <= loglikstore) # 0.6324

hist(loglikstore)
abline(v=loglik, col="red", lwd=2)

pdf(paste("balance_interactions.pdf"),w=6,h=5)
hist(loglikstore, breaks = 100, main = paste("Sampling distribution of simulated log likelihoods"), xlab=("Log Likelihood"))
abline(v=loglik, col="red", lwd=3)
dev.off()